
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: aarch64-apple-darwin20 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> # Script to run QAP
> library(sna)
Loading required package: statnet.common

Attaching package: ‘statnet.common’

The following objects are masked from ‘package:base’:

    attr, order

Loading required package: network

‘network’ 1.18.1 (2023-01-24), part of the Statnet Project
* ‘news(package="network")’ for changes since last version
* ‘citation("network")’ for citation information
* ‘https://statnet.org’ for help, support, and other information

sna: Tools for Social Network Analysis
Version 2.7-1 created on 2023-01-24.
copyright (c) 2005, Carter T. Butts, University of California-Irvine
 For citation information, type citation("sna").
 Type help(package="sna") to get started.

> library(doParallel)
Loading required package: foreach
Loading required package: iterators
Loading required package: parallel
> library(doRNG)
Loading required package: rngtools
> 
> # Load the followers adjacency matrix
> y <- readRDS("processed_data/mentions_adjacencyMatrix.Rds")
> # Load the predicting matrix 
> load("processed_data/QAP_predicting_matrices.RData")
> 
> diag(y) <- NA
> 
> y_v <- c(y)
> 
> 
> Var.names <-c("State_Similarity",
+               "Party_Similarity",
+               "Chamber_Similarity",
+               "Gender_Similarity",
+               "Race_Similarity",
+               "Difference_in_Legislatures_Profeshionalism",
+               "Dem_Sender_Effect",
+               "Rep_Sender_Effect",
+               "House_Sender_Effect",
+               "Female_Sender_Effect",
+               "Profesh_Sender_Effect",
+               "Black_Sender_Effect",
+               "Latino_Sender_Effect",
+               "Asian_Sender_Effect",
+               "Mena_Sender_Effect",
+               "Multi_Sender_Effect",
+               "Native_Sender_Effect",
+               "Democrat_Receiver_Effect",
+               "Republican_Receiver_Effect",
+               "House_Receiver_Effect",
+               "Female_Receiver_Effect",
+               "Profesh_Receiver_Effect",
+               "Black_Receiver_Effect",
+               "Latino_Receiver_Effect",
+               "Asian_Receiver_Effect",
+               "Mena_Receiver_Effect",
+               "Multi_Receiver_Effect",
+               "Native_Receiver_Effect",
+               "Same_Party_Same_State",
+               "Same_Chamber_Same_State",
+               "Same_Gender_Same_State",
+               "Same_Race_Same_State",
+               "Contiguous_States")
> 
> for(i in 1:length(Var.names)){
+   xi <- predicting_matrices[i,,]
+   diag(xi) <- NA
+   print(sum(is.na(xi)))
+   if(i == 1){
+     xdf <- cbind(c(xi))
+   }
+   if(i > 1){
+     xdf <- cbind(xdf,cbind(c(xi)))
+   }
+ }
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
[1] 4108
> 
> send_id <- matrix(1:nrow(y),nrow(y),nrow(y))
> rec_id <- t(send_id)
> diag(send_id) <- NA
> diag(rec_id) <- NA
> 
> xdf <- cbind(xdf,c(send_id),c(rec_id))
> 
> xdf <- data.frame(xdf)
> 
> names(xdf) <- c(Var.names,c("send_id","rec_id"))
> 
> xydf <- data.frame(y_v,xdf)
> 
> xynd <- na.omit(xydf)
> 
> covs <- paste(Var.names,collapse="+")
> 
> names(xynd)[1] <- "yij"
> 
> library(fastglm)
Loading required package: bigmemory
> system.time(est0 <- fastglm(as.matrix(xynd[,Var.names]),xynd$yij,family=gaussian(),method=3))
   user  system elapsed 
  6.163  12.746  33.464 
> 
> 
> load("simulated_mention_networks.RData")
> rm(xdf)
> gc()
             used    (Mb) gc trigger    (Mb) limit (Mb)   max used    (Mb)
Ncells     665724    35.6    1236465    66.1         NA    1236465    66.1
Vcells 2777423045 21190.1 3597346796 27445.6      32768 2996423145 22860.9
> rm(predicting_matrices)
> gc()
             used    (Mb) gc trigger    (Mb) limit (Mb)   max used    (Mb)
Ncells     665717    35.6    1236465    66.1         NA    1236465    66.1
Vcells 2220526126 16941.3 3597346796 27445.6      32768 2996423145 22860.9
> rm(xynd)
> gc()
             used    (Mb) gc trigger    (Mb) limit (Mb)   max used    (Mb)
Ncells     665667    35.6    1236465    66.1         NA    1236465    66.1
Vcells 1621581707 12371.7 3597346796 27445.6      32768 2996423145 22860.9
> 
> library(biglm)
Loading required package: DBI
> 
> estim <- NULL
> twop <- NULL
> for(i in 1:50){
+   
+   # qap unpermuted
+   ys <- sim_amats[[i]]
+   
+   diag(ys) <- NA
+   
+   y_vs <- c(y)
+   
+   xydf$y_vs <- y_vs
+   
+   covs <- paste(Var.names,collapse="+")
+   
+   form <- as.formula(paste("y_vs~",covs,sep=""))
+   
+   system.time(est_true <-  biglm(form,data=xydf))
+   
+   filenm =paste("./qap_power_results/qap_power_results/qap_sim_mention",i,".RData",sep="")
+   load(filenm)
+   
+   obs_ts <- summary(est_true)$mat[,1]/summary(est_true)$mat[,4]
+   
+   pgabs <- NULL
+   for(j in 1:length(coef(est_true))){
+     pgabs <- c(pgabs,mean(abs(perm_coef[,j]) > abs(obs_ts[j])))
+   }
+   
+   estim <- cbind(estim,coef(est_true))
+   twop <- cbind(twop,pgabs)
+   
+   print(i)
+ }
[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
[1] 10
[1] 11
[1] 12
[1] 13
[1] 14
[1] 15
[1] 16
[1] 17
[1] 18
[1] 19
[1] 20
[1] 21
[1] 22
[1] 23
[1] 24
[1] 25
[1] 26
[1] 27
[1] 28
[1] 29
[1] 30
[1] 31
[1] 32
[1] 33
[1] 34
[1] 35
[1] 36
[1] 37
[1] 38
[1] 39
[1] 40
[1] 41
[1] 42
[1] 43
[1] 44
[1] 45
[1] 46
[1] 47
[1] 48
[1] 49
[1] 50
> 
> psig <- NULL
> for(i in 1:length(coef(est0))){
+   same_sign <- sign(estim[i+1,]) == sign(coef(est0))[i]
+   sig <- mean((twop[i+1,] <= 0.05)*same_sign)
+   psig <- c(psig,sig)
+ }
> 
> # calculate standardized coefficients
> abs_std_coef <- coef(est0)/est0$se
> abs_std_coef <- abs(abs_std_coef)
> 
> names(coef(est0))[abs_std_coef > 2 & psig < 0.8]
[1] "Dem_Sender_Effect"        "Profesh_Sender_Effect"   
[3] "Democrat_Receiver_Effect" "Profesh_Receiver_Effect" 
> 
> library(ggplot2)
> # Basic scatter plot
> df <- data.frame(power = psig,zstat = abs_std_coef,
+                  col = ifelse(abs_std_coef > 2 & psig < 0.8,"red","blue"),
+                  shp = ifelse(abs_std_coef > 2 & psig < 0.8,17,19))
> 
> ggplot(df, aes(x=zstat, y=power)) +
+   geom_point(size=2,shape =df$shp,colour=df$col) + 
+   scale_x_continuous(trans='log2') +
+   geom_hline(yintercept=0.8) +
+   geom_vline(xintercept=2) +
+   theme_bw()
> 
> ggsave("power_results_mention.pdf",width=4,height=3.5)
> underpowered <- names(coef(est0))[abs_std_coef > 2 & psig < 0.8]
> write.csv(underpowered,file="underpowered_mention.csv")
> 
> 
> 
> 
> 
> 
> 
> proc.time()
    user   system  elapsed 
 696.964  581.686 3354.125 
